Funnel landscape and mutational robustness as a result of evolution under thermal 
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Using a statistical-mechanical model of spins, evolution of phenotype dynamics is studied. Con- 
figurations of spins and their interaction J represent phenotype and genotype, respectively. The 
fitness for selection of J is given by the equilibrium spin configurations determined by a Hamiltonian 
with J under thermal noise. The genotype J evolves through mutational changes under selection 
pressure to raise its fitness value. From Monte Carlo simulations we have found that the frustration 
around the target spins disappears for J evolved under temperature beyond a certain threshold. 
The evolved Js give the funnel-like dynamics, which is robust to noise and also to mutation. 

PACS numbers: 87.10.-c,87.10.Hk,87.10.Mn,87.10.Rt 



Under fixed conditions, biological systems evolve to 
increase their fitness, determined by a biological state 
-phenotype- that is shaped by a dynamical process. This 
dynamics is generally stochastic as they are subject to 
thermal noise, and the rule for the dynamics is controlled 
by a gene that mutates through generations. Those genes 
that produce higher fitness values have a higher chance 
of survival. For example, the folding dynamics of a pro- 
tein or t-RNA shapes a structure under thermal noise 
to produce a biological function, while the rule for the 
dynamics is coded by a sequence of DNA. 

The phenotype that provides fitness is expected to be 
rather insensitive to the type of noise encountered dur- 
ing the dynamical process [l|, l3, ISl , to continue producing 
fitted phenotypes. The first issue is to determine what 
type of dynamics is shaped through evolution to achieve 
robustness to noise. To have such robustness it is prefer- 
able to utilize a dynamical process that produces and 
maintains target phenotypes of high fitness smoothly and 
globally from a variety of initial configurations. In fact, 
the existence of such a global attraction in theprotein 
folding was proposed as a consistency principle [4| and a 
funnel landscape [5[, while similar global attraction has 
been discovered in gene regulatory networks [6|. Despite 
the ubiquity of such funnel landscapes for phenotype dy- 
namics, little is understood on how these structures are 
shaped by the evolution [Tl, la, lof . The second issue we ad- 
dress is the conditions under which funnel-like dynamics 
evolves. 

Besides its robustness to noise, the phenotype is also 
expected to be robust against mutations in the genetic 
sequence encountered through evolution. Despite recent 
studies suggesting a relationship between robustness to 
thermal noise and robustness to mutation [l|, lid . I 111 . Il2l . 
13| , a theoretical understanding of the evolution of the 
two is still insufficient. The question of whether robust- 
ness to thermal noise also leads to robustness to mutation 
constitutes the third issue discussed in this paper. All 
these questions concerning robustness to noise and muta- 



tion or the shaping of the funnel-like landscape can be an- 
swered by introducing an abstract statistical-mechanical 
model of interacting spins, whose Hamiltonian evolves 
over generations to achieve a higher fitness. 

Let us consider a system of N Ising spins interacting 
globally. In this model, configurations of spin variables Si 
and an interaction matrix Jy with i, j ~ I, ■ ■ ■ , N repre- 
sent phenotype and genotype, respectively. For simplic- 
ity, Jij is restricted to the two values ±1 and is assumed 
to be symmetric: Jij ~ Jji. The dynamics of the phe- 
notype denoted by S is given by a flip-flop update of 
each spin with an energy function, defined by the Hamil- 
tonian H{S\J) = ^~m J2i<j JijSiSj, for a given set of 
genotypes denoted by J. We adopt the Glauber dynam- 
ics as an update rule, where the N spins are in contact 
with a heat bath of temperature Ts- After the relaxation 
process, this dynamics yields an equilibrium distribution 



for a given J, P{S\J,Ts) = e 



-l3sH{S\J) 



/ZsiTs), where 



(3s = l/Ts and Zs{Ts) = TVsexp[-/3sF(5| J)]. 

The genotype J is transmitted to the next generation 
with some variation, whereas genotypes that produce a 
phenotype with a higher level of fitness are selected. We 
assume that fitness is a function of the configuration of 
target spins t, a given subset of S, with size t. The time- 
scale for genotypic change is generally much larger than 
that for the phenotypic dynamics, so that the variables 
S are well equilibrated within the unit time-scale of the 
slower variable J. Such being the case, fitness should be 
expressed as a function of the target spins S averaged 
with respect to the distribution. Considering the gauge 
transformation llSl. we define fitness as 



*(J|r5) 



n ^(^^ 
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Sj) 



without losing generality, where (■ • ■ ) denotes the expec- 
tation value with respect to the equilibrium probabil- 
ity distribution. In other words, the expected fitness is 
given by the probability of the "target configuration" in 




FIG. 1: (color online) Contour maps on Ts and Tj of (a) 
fitness, (b) energy, and (c) $2, where 1 — $2 is the frustration 
around the target spins, for evolved J at a given Ts and Tj 
(see text for details). A^ = 15 and t — 3. For each generation 
of the genotype dynamics, the average in equilibrium is taken 
over 1500 MC steps after discarding the first 1500 MC steps, 
which are sufficient for equilibration. The data are averaged 
over the last 10'' generations. 



which all target spins are aligned parallel under equilib- 
rium conditions. Note that in our model, only the target 
spins contribute to the fitness, and as a result, the spin 
configuration for a given fitness value has redundancy. 

The genotype J evolves as a result of mutations, ran- 
dom flip-flop of the matrix element, and the process of 
selection according to the fitness function. We again 
adopt Glauber dynamics by using fitness instead of the 
Hamiltonian in the phenotype dynamics, where J is in 
contact with a heat bath whose temperature Tj is dif- 
ferent from Ts- In particular, the dynamics is given 
by a stochastic Markov process with the stationary dis- 
tribution P{J^Ts,Tj) ^ eP-''^'^-^'^s)iZj{Ts,Tj), where 
pj = 1/Tj and Zj{Ts,Tj) = Trjcxp[A7*(J, Ts)]- Ac- 
cording to the dynamics, genotypes arc selected some- 
what uniformly at high temperatures Tj, whereas at low 
Tj^ genotypes with higher fitness values are preferred. 
The temperature Tj represents selection pressure. 

Next, we study the dependence of the fitness and en- 
ergy on Tg and Tj, given by ^(Ts, Tj) = {■^{J\Ts)).j and 
E{Ts,Tj) ^ {{HiS\J)))j, respectively, where (• • ■ )j de- 
notes the average with respect to the equilibrium proba- 
bility distribution, P{J,Ts,Tj). For the spin dynamics 
(unless otherwise mentioned), the exchange Monte Carlo 
(EMC) simulation [ij] is used to accelerate the relaxation 
to equilibrium. Indeed, we have confirmed the equilib- 
rium distribution for the simulations below. Two pro- 
cesses are carried out alternately: the equilibration of S 
with the EMC and the stochastic selection of J according 
to the fitness value estimated through the first process. 

Fig. [1] (a) and (b) show dependence of the fitness and 
energy on Ts and Tj, respectively, for TV = 15 and i = 3. 
For any Ts, the fitness value decreases monotonically 
with Tj. However, Ts influences the slope of the de- 
crease significantly. The fitness for sufficiently low Ts 
remains at a high level and decreases only slightly with 
an increase in Tj, while for a medium value of Ts, the 
fitness gradually falls to a lower level as a function of 
Tj and eventually, for a sufficiently high value of Ts, it 



never attains a high level. This implies that the struc- 
ture of the fitness landscape depends on Ts, at which the 
system has evolved. The energy function, on the other 
hand, shows a significant dependence on Ts- Although 
the energy increases monotonically with Ts for high Tj, 
it exhibits non-monotonic behavior at a low Tj and takes 
a minimum at an intermediate Ts- The J configurations 
giving rise to the highest fitness value generally have a 
large redundancy. At around Ts — 2.0, using a fiuc- 
tuation induced by Ts, a specific subset of adapted J's 
giving rise to lower energy is selected from the redundant 
configurations with higher fitness. 

In the medium-tcmpcrature range, such Js that yield 
both lower energy and higher fitness are evolved. Now 
we study the characteristics of such Js. According to 
statistical physics of spin systems, triplets of interactions 
that satisfy JijJjkJki < are known to yield frustra- 
tion, which is an obstacle to attaining the unique global 
energy minimum [l5|. In our model, however, the tar- 
get spins play a distinct role. Hence, it becomes nec- 
essary to quantify the frustration by distinguishing tar- 
get and non-target spins. In accordance with the "fer- 
romagnetic" fitness condition for target spins, we define 
$1 as the frequency of positive coupling among target 

spins, ^i{Ts,Tj) = ^^^(j2,^,^tJ^])y Under ferro- 
magnetic coupling, the target configurations are energet- 
ically favored, i.e., $1 = 1, for which no frustration exists 
among the target spins. Next, for a measure of the de- 
grees of frustration between target and non-target spins, 
and among non-target spins themselves, we define $2 and 
$3 as 



^2{Ts,Tj 
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where Cj ~ is the total number of possible pairs among 
the non-target spins. Here, $2 is the fraction of the in- 
teraction pairs between target and non-target spins that 
satisfy JikJkj = 1 (i, j € *, k ^ t). If $2 = 1, no frus- 
tration is introduced by the interaction between target 
and non-target spins, so that the energy minimum of the 
target configuration is conserved by such interaction. If 
<1>3 = 1, the target configurations do not introduce frus- 
tration in Jkis {k, I ^ t). If $1 = $2 = <1>3 = 1, there is 
no frustration at all over the interactions, as suggested 
by the Mattis model [17| , which can be transformed into 
ferromagnetic interactions by gauge transformation [l5|. 
For J that is evolved under given Ts and Tj, we have 
computed $i,$2 and <1>3. Fig. [2] shows the dependence 
of $1, $2 and $3 on Ts at a fixed T^ = 0.5 x IQ-^. For 



Ts > Tf-, (f>i takes the value -^1 [l6|, so that a target 
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FIG. 2: (color online) Dependence of local frustrations on 
Ts,; $i(n),<l?2(0)5'3E>3(A) represents the left axis and /ic(*) 
represents the right axis by fixing Tj at 0.5 x 10~ . The 
data are computed by taking averages over 150 genotypes J 
evolved at a given temperature Ts. /Jc is a threshold value, 
beyond which the fitness of the mutated genotype begins to 
decrease (see Fig. |3|. The transition points T§^ and Tg^ are 
estimated as a temperature at which $2 deviates from 1. 



FIG. 3: (color online) Relaxation dynamics of the averaged 
magnetization of target spins, {mt)o on the adapted interac- 
tion J^f' for Ts = 10"^ and Ts = 2.0. The magnetization 
{mt)o is evaluated by taking the average over 30 initial con- 
ditions for each J^ ^ and 1000 different samples of J^ ''. The 
inset shows the dependence of the estimated convergent value 
of {mt)o on Ts, m* represents the right axis and the relax- 
ation time T represents the left axis. 



configuration is embedded as an energetically favorable 
state, while no specific patterns, apart from the target, 
are embedded in the spin configuration 



where (• • ■ )o denotes the average over the randomly cho- 
sen initial conditions. As shown in Fig. [31 the relaxation 
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For rgi <Ts< Tf, $2 also takes the value around 1, P™*^^^^ ^""^ ^Ts' evolved at low temperatures is much 



implying that frustration is not introduced by means of 
interactions with a non-target spin. In this temperature 
range, $3 is not equal to 1, except for Ts ~ 2.0 where 
the Mattis state is shaped. When $2 ~ 1 and $3 ^ 1, 
frustration is not completely eliminated from the non- 
target spin interactions, in contrast to the Mattis state. 
Here, such a J configuration without frustration around 
the target spins (but with frustration between non-target 
spins) is referred to as "local Mattis state" (LMS), as 
characterized by <i>i = (f>2 "^ 1 and <I>3 7^ 1. The inter- 
actions J that form such LMSs arise as a result of the 
evolution at T|^ < Tg < T|^, where both a fitted target 
configuration and a lower energy level are achieved. As 
shown in Fig. [TJc), the Ts range in which LMS is shaped 
becomes narrower with an increase in Trjl8|. For suffi- 
cient low Tj, there are three phases: Ts < Ts^(Tj), the 
phase in which frustration remains in spite of adapta- 
tion; T^^{Tj) <Ts< Tf{Tj), the phase giving LMSs; 
Ts > Tg^{Tj), the phase in which no adaptation and 
frustration is seen. 

Let us consider the relaxation dynamics of spins for 
each J adapted through evolution under a given Ts and 
Tj, denoted as J^g^'^ where Tj is fixed at 0.5 x 10"^. 
To understand how the relaxation dynamics depends on 
J^s^, instead of using EMC, we adopt standard MC by 
with the temperature Tg, fixed at 10~^, independently 

''° ^ We compute the temporal 
^_^^^S,\. Re- 
laxation dynamics of {mt)o for J°^p, Ts = 10-3(< T^i), 
and Ts = 2.0{T§^ < Ts < Tf) are plotted in Fig. El 



of Ts used in obtaining j!^^' 

change of the target magnetization mt = | X^iet ' 



slower. Furthermore, {mt)o converges to a value ttij lower 
than 1 and remains at that value for a long time. Depend- 
ing on the initial condition, the spins are often trapped 
at a local minimum, so that the target configuration is 
not realized over a long time span. 

Such dependence on initial conditions is not observed 
for J^ ^ for Ts > Ts^, where {mt)o approaches 1 some- 
what quickly. From an estimate of the convergent value 
of the target magnetization ml within the above MC 
time scale, we obtain the relaxation time r by fitting to 
the function (mt)o(s) = ml -|-cexp(— s/r), where s is the 
MC step of the spin dynamics. The parameters m.^ and t 
are plotted against Ts in the inset of Fig. [3l which shows 
the increase of r and the decrease of m.^ from 1 with the 
decrease of Ts below Tg^ . These results imply that the 
energy landscape for the interaction J^ ^ is rugged for 
Ts < Tg^, as in a spin-glass phase, whereas it is smooth 
around the target for T^^ < Ts < Tf. Thus, this land- 
scape is interpreted as a typical funnel landscape. It 
demonstrates a transition from the spin-glass phase to 
the funnel at Tf (see also 01). 

Now, in the evolved genotypes, let us examine the ro- 
bustness that represents the stability of J's fitness with 
respect to changes in the J configuration. From the 
adopted genotype J^ ^ , mutations are imposed by fiip- 
flopping the sign of a certain fraction of randomly cho- 
sen matrix elements in J^ ^ . The value of the fraction 
represents the mutation rate ^. We evaluate the fit- 



ness of the nrutated Jt^^(m) at T'g = 10'^ {^ Ts), i.e., 
^(J^ ^{^J)\Ts = 10~^), by taking an average over 150 
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FIG. 4: (color online) Fitness of the mutated J as a function 
of the mutation rate ^ for Tg — 10"'* (solid curve) and Ts ~ 
2.0 (dotted curve). For each adapted genotype, the mutated 
genotypes J are generated 150 times by flipping randomly 
chosen elements Jij with J governed by the rate /i, and the 
average is taken over 150 adapted genotypes. 



samples of mutated J^ ^(^f)- Fig. [5] shows /i dependence 
of the fitness for Ts = 10^^ and Ts = 2.0. For low values 
of T5, the fitness of mutated J^ ^{lA exhibits a rapid 
decrease with respect to the mutation rate, but for Ts 
between Tg^ and T^^, it does not decrease until the mu- 
tation rate reaches a specific value. We define ndTs) as 
the threshold point in the mutation rate beyond which 
the fitness ^(J^ ^{fJ-)\T'g) begins to decrease from unity. 
Fig. [21 shows the dependence of ^c on Ts, which has a 
plateau at T^^ < Ts < Tf where $2 is unity This 
range of temperatures that exhibits mutational robust- 
ness agrees with the range that gives rise to the LMS. In 
other words, mutational robustness is realized for a geno- 
type with no frustration around the target spins. Evo- 
lution in a mutationally robust genotype J is possible 
only when the phenotype dynamics is subjected to noise 



within the range Tg 



"^ <Ts < Tf. This mutational ro- 



bustness is interpreted as a consequence that the fitness 
landscape becomes non- neutral for Ts > T|^ [18j. 

To check the generality of the transition to the LMS 
as well as the mutational robustness, we have examined 
the model by increasing the number of target spins i, and 
confirmed that the LMSs evolve at an intermediate range 
of Ts (that depends on t) , where both lower energy and 
higher fitness are realized together with mutational ro- 
bustness, while the actual fitness value therein decreases 
with an increase in t. Simulations with larger N (up to 
30) have also confirmed the evolution of the LMSs at an 
intermediate range of Ts [18[ . 

In this study, in order to elucidate the evolutionary 
origin of robustness and funnel landscape, we have con- 
sidered the evolution of a Hamiltonian system to gener- 
ate a specific configuration for target spins. The findings 
can be summarized as follows. First, as a result of the 
formation of a funnel landscape through the evolution 
of the Hamiltonian. robustness to guard against noise is 



achieved in the dynamic process. Such shaping of dynam- 
ics is possible only under a certain level of thermal noise, 
given by temperatures T|^ < Ts < T|^. Second, under 
such a temperature range in the process, a funnel-like 
landscape that gives rise to a smooth relaxation dynam- 
ics toward the target phenotype is evolved to avoid the 
spin-glass phase. This may explain the ubiquity of such 
funnel-type dynamics observed in evolved biological sys - 
tems such as protein folding and gene expression p, \l^ . 
Third, this robustness to thermal noise induces robust- 
ness to mutation; this observation has also been discussed 
for gene transcription network models [I4I . Relevance of 
thermal noise to robust evolutions is thus demonstrated. 

The funnel-like landscape evolved at Tji < Ts < Tf 
is characterized by the local Mattis state without frustra- 
tion around the target spins. This allows for a smooth 
and quick relaxation to the target configuration, which is 
in contrast with the relaxation on a rugged landscape in 
spin-glass evolved at Ts < Tg^ , where relaxation is often 
trapped into metastable states. Theoretical analysis of 
random spin systems such as replica symmetry breaking 
will be relevant to the local Mattis state and mutational 
robustness [18[ , as the model discussed in this letter is a 
variant of spin systems with two temperatures [l9| . 
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